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ABSTRACT 


QR  is  the  standard  method  for  finding  all  the  eigenvalues  of  a  symmetric  tridiagonal  matrix.  It  produces  a 
sequence  of  similar  tridiagonals.  It  is  well  known  that  the  QR  tran^ormation  from  T  to  f  is  backward 
stable.  That  means  that  the  computed  f  is  exactly  orthogonally  similar  to  a  matrix  cbse  to  T.  It  is  also 
known  that  the  algorithm  sometimes  exhibits  forward  instability.  That  means  that  the  computed  f  is  not 
close  to  the  exact  f. 

For  the  purpose  of  computing  eigenvalues  the  property  of  backward  stability  is  all  that  one  requires.  How¬ 
ever  the  QR  transformation  has  other  uses  and  there  forward  stability  is  wanted. 

This  report  analyzes  the  forward  instability  and  shows  that  it  occurs  only  when  the  shft  causes  premature 
deflation.  We  show  that  forward  stability  is  governed  by  the  behavior,  in  exact  arithmetic,  of  a  pair  of  vari¬ 
ables  and  we  establish  tight  upper  and  lower  bounds  on  their  derivatives  with  respect  to  change  in  the  shift 
parameter. 
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/  1.  Summary  and  Notation 


1 

i. 


Ll  Introductioii 

.This  study  is  in  the  area  of  matrix  eigenvalue  computations. 


^  im: 

<^Thc 


lie  Householder-QR  algorithm  has  become  the  standard  method  for  diagonalizing  a  symmetric 
matrix.  First  the  matrix  is  reduced  to  tridiagonal  form  7  by  a  technique  introduced  by  A.  Householder  in 
1958.  Next  the  tridiagonal  matrix  T  is  ^gontdized  byLSUccessiye.  applications  of  the  QR  transfcxmation 
with  shifts.  Moreover  it  is  well  known4<-ae«  f  Wilkinson,  chap:-  3-]  >that  the  QR  transformation  is  back¬ 
ward  stable.  That  means  that  the  computed  transform  is  exactly  orthogonally  similar  to  a  matrix  close  to  the 
old  one.  It  is  also  known  to  the  expert^C~5ec  tWllkinson,  4958  },  [  Golub  &  Kahan  ]  ^  [  Stewart,  1970  ] ) 
-that  the  transformation  sometim^  exhibits  forward  instability.  That  means  that  the  computed  output  is 
far  firtxn  the  result  obtained  with  exact  arithmetic.  For  the  purpose  of  computing  eigenvalues  and  eigenvec¬ 
tors  the  property  of  backward  stability  is  all  that  one  requires.  However  the  QR  transformation  has  other 
uses  and  in  some  cases  forward  stability  is  desirable. 

rqx>rt  analyzes  the  forward  instability  of  QR  and  shows  that  it  occurs  only  when  the  shift  is 
very  close  to  eigenvalues  with  a  special  property.^or  some  matrices  there  may  be  no  eigenvalues  with  this 
property  and  in  such  cases  the  algorithm  is  forward  stable  for  any  value  of  the  shift.  The  study  was 
prompt^  by  the  discovery  that  the  dominant  facto^  determining  inst^ility  is  the  sensitivity  of  several  key 
quantities  to  small  changes  in  the  shift  parameter,  '^s  sensitivity  dominates  the  effect  of  perturbation  in  ail 
the  other  variables.  i  — ^  ^ 

The  technical  contributions  of  this  study  are;  \  y  - 

1)  The  observation  that  instability  is  equivalent  to  premature  deflation  that  occurs  when  the 
shift  is  almost  an  eigenvalue  of  several  consecutive  submatrices; 

2)  Establishing  the  central  role  of  a  certain  pair  of  variables  associated  with  each  plane  rota¬ 
tion  used  in  the  algorithm.  It  is  the  evolution  of  the  norm  of  this  2-vector  that  governs  the 
accuracy  of  the  computed  angles. 

3)  Bounds  on  the  derivatives  of  the  cosines  of  the  rotation  angles  and  other  key  quantities. 

4)  Upper  and  lower  bounds  on  the  last  components  of  normalized  eigenvectors  in  terms  of 
eigenvalues. 

1.2  Organization  of  this  study 

This  first  secticm  gives  the  summary  of  this  study.  Along  with  it,  the  notation  conventions  to  be  used^ 
in  the  discussions  is  presented. 

In  section  2,  the  07  transformation  is  defined  and  several  examples  are  exhibited  to  show  that  some, 
times  the  QR  transformation  ot  T  is  forward  stable  and  sometimes  it  is  not  Also  in  section  2,  severaK 
needed  spectral  properties  of  T  are  described.  An  useful  result  is  Theorem  23  that  gives  upper  and  lower 
bounds  on  the  last  element  of  a  normalized  eigenvector  of  T. 

In  section  3,  the  implementation  of  the  QR  transformation  on  T  is  discussed.  In  the  process,  the 
important  intermediate  quantities  are  introduced  along  with  the  relations  among  them. 


For 

In  section  4,  these  quantities  are  analyzed  in  terms  of  the  shift  a.  The  main  results  of  this  effort  are  — 
1)  the  derivatives  of  ttt ,  c/t+i  ( to  be  defined  in  (3.1.2)  )  will  attain  their  extrema  at  the  eigenvalues  of  ;  2)  ^ 
if  some  eigenvalue  of  T*  is  very  close  to  some  eigenvalue  of  Tk-u  the  absolute  values  of 
ditkida,  dckida  for  a  in  the  vicinity  of  this  special  X^*^  can  be  large. 

In  section  5,  the  intermediate  quantities  produced  in  finite  precision  arithmetic  by  a  particular  QR 
implementation  are  formulated  into  a  convenient  matrix-vector  form  and  the  influence  of  the  roundoff 
errors  is  represented  by  a  tridiagonal  error  matrix.  The  important  usage  of  this  error  matrix  is  in  the  Per¬ 
turbed  Commutative  Law  which  is  used  to  explain  the  forward  instability  of  TQR . 
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U  Notation 

Throughout  this  woilc  we  will  use  T  to  denote  the  real  n  x  n  symmetric  tridiagonal  matrix.  In  general, 
upper  case  Roman  letters  will  denote  matrices,  lower  case  Bold  letters  will  denote  vectors  and  lower  case 
Roman  letters  will  denote  scalars.  Upper  case  Greek  letters  are  used  for  special  matrices  (  usually  diagonal 
).  lower  case  Greek  letters  are  also  scalars. 

The  matrix  denotes  the  kxk  identity  matrix.  The  matrix  denotes  the  k-th  leading  principal 
submatrix  of  T.  The  norm  II  ■  II  is  the  Euclidean  norm.  The  tran^se  operation  is  denoted  by '  ( e.g.  M'  is 
read  M  transpose  ). 

1.4  ( More  technical )  Summary 

The  actual  algorithm  used  to  implement  the  QR  transformation  employs  a  sequence  of  plane  rota¬ 
tions  in  planes  (1,2),  (2, 3),  (3,4), (n-l,n  )  that  change  T,  into  T,  by  chasing  a  coiain  l^ge  in  the 
tridiagonal  fonn  down  the  matrix  and  off  the  end.  The  proof  that  this  process  is  equivalent  to  the  formal 
definition  of  the  QR  transform  (  namely  (i)T  -al  =  QR;  (ii)  f  =RQ  +a/, )  dq)ends  strongly  on  the 
unreduced  property  ( to  be  defined  in  section  2.2  )  of  7.  In  a  finite  precision  environment  we  must  antici¬ 
pate  a  divergence  of  these  two  processes  when  any  subdiagonal  elements  |3t(2^jk^n)are  small 
enough,  except  for  the  last  one.  In  other  words  it  is  not  surprising  that  small  ^k>(.k  <n  )  causes  fraward 
instability  in  QR . 

Nevertheless  this  possible  closeness  to  leducibility  is  not  the  only  cause  of  instability.  The  purpose  of 
this  study  is  to  elucidate  exactly  how  forward  instability  can  occur  both  with  and  without  small  pt, 
(k  <n). 

Forward  instability,  if  it  occurs,  is  quite  dramatic.  The  cosines  of  the  rotation  angles,  cj,  C2 . c,-i 

exhibit  the  following  t^avior.  Up  to  some  k<n-l  the  computed  and  the  exact  c,-  have  the  same 
exponents  which  eventually  diminish  to  log  {roundoff  unit)H.  For  itk  the  true  I  Ci  I  continues  to 
decrease  while  the  comput^  I  Cj  I  increases  in  such  a  way  that  the  products  I  Cj  ^  1=0 (e).  This  holds 
until  I  Cj  I  <  e.  As  I  Cj  I  diminishes  even  further,  I  Ci  I  stays  close  to  1  unless  Pi  suddenly  drops.  An  iso¬ 
lated  vanishing  of  Ci  (i  <k)  does  no  harm. 

One  result  of  our  study  is  that  forward  instability  is  always  associated  with  "premature  deflation".  In 
the  scenario  given  in  the  previous  paragraph  it  h^^ns  that  after  rotation  k  the  elements  {k-l,k), 
{k,  k-l ),  (i ,  k+1 ),  {k+l.  A; )  are  all  on  the  order  of  II  T  II .  The  shift  appears  in  position  {k,k)  and  is 
correct  to  workitjg  precision.  If  row  k  and  column  k  are  deleted  from  the  new  matrix  to  obtain  then  the 
eigenvalues  of  7^*^  will  be  the  remaining  eigenvalues  of  7.  to  within  working  accuracy.  If  the  shift  is  a 
well  isolated  eigenvalue  of  7^  then  its  eigenvector  can  be  constructed  from  the  rotation  angles  up  to  . 

The  occurrence  of  forward  instability  is  not  connected  with  the  presence  of  clusters  of  close- 
eigenvalues  in  7,.  It  is  caused  by  the  shift  being  an  eigenvalue  of  and  Tt  (to  working  accuracy ).  It  so 
happens  that  when  this  occurs  the  shift  will  be  an  eigenvalue  of  all  the  principal  submatrices  7t_i,  7*,  7t+i, 
...,  7, ,  ( to  working  accuracy ). 

We  have  described  instability  in  terms  of  the  cosines  c)  because  they  are  more  familiar.  However  a 
better  indicator  of  stability  is  (  where  n,  is  one  of  the  variables  that  appears  in  TQR .  See 

section  5  for  more  details. 

1.5  Application  to  the  Lanezos  Algorithm 

In  general  the  algorithm  TQR  ( tridiagonal  QR  )  is  forward  stable  for  all  choice  of  shifts.  However 
there  is  an  important  triplication  where  instability  is  endemic  and  it  was  the  gradual  realization  of  this 
uncomfortable  fact  that  1^  to  our  study. 

The  Lanezos  itoation  produces  a  symmetric  tridiagonal  matrix  to  which  a  new  row  and  column  are 
added  at  each  step.  At  the  end  of  s/ep  k  the  algorithm  has  produced  tridiagonal  7^  and  an  extra  number 
P*+i-  7*  represents  the  projection  of  some  given  linear  operator  A  on  a  special  -dimensional  subspace. 
These  growing  tridiagonals  are  special  because,  as  k  increases,  some  eigenvalues  stabilize.  In  other  words 
an  eigenvalue  of  Tt  appears  to  be  equal  ( to  working  i^-ecisiem  )  to  an  eigenvalue  of  Tt+i  and  to  an  eigen¬ 
value  of  Tt^.2.  and  so  on.  These  stabilized  values  are  eigenvalues  of  the  linear  operator  A  on  R".  In  an 
implementation  of  the  Lanezos  algorithm  it  is  convenient  to  get  rid  of  these  converged  eigenvalues  by 
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defladng  them  from  T,.  At  first  sight  this  ^)pears  to  be  impossible  because  at  step  k  (k  <n  )\he  matrix 
r,  is  not  fully  known.  However  this  deflation  is  possible  even  occurs  naturally  ( in  exact  arithmetic  ) 
when  the  QR  algorithm  is  tqiplied  to  Th+i  with  the  correct  shift  The  unknown  component  in  position 
(k+1,  k+l )  is  not  altered  when  the  QR  transformation  is  halted  at  step  k.  It  is  only  necessary  to  delete 
row  and  column  k . 

What  happened  to  us  was  that  we  did  not  always  know  the  ctxrect  value  k .  When  the  QR  transfor¬ 
mation  was  forced  to  continue  beyond  the  right  place  then  the  results  were  terrible.  As  a  result,  the 
deflation  by  QR  failed  and  the  resulting  tridiagonals  were  wrong.  Of  course  deflation  had  occurred  earlier 
but  we  did  not  look  for  it  The  process  had  encounted  forward  instability.  In  the  spirit  of  knowing  your 
enemy  this  investigation  was  launched. 


2.  Background  Information 

This  section  covers  the  definition  of  the  QR  transformation  and  its  relation  to  eigenvalue  deflation 
and  eigenvector  calculation  for  a  symmetric  tridiagonal  matrix  T.  The  examples  in  section  23  show  that 
the  QR  transformation  on  T  can  sometimes  be  violently  unstable  in  the  forward  sense.  For  easy  reference, 
several  needed  spectral  properties  of  T  are  collected  into  section  2.4  from  [  Parlett,  chap.7  ].  Theorem  2.3 
gives  upper  and  lower  bounds  on  the  bottom  element  of  a  nonnalized  eigenvector. 

2.1  The  QR  transfwmation 

This  well  established  procedure  is  described  in  several  books;  e.g.  [  Wilkinson,  chap.8  ],  [  Stewart, 
chap.7  ],  [  Parlett,  chap.8  ],  [  Golub  &  Van  Loan,  chap.7  ].  Here  we  will  reproduce  only  what  we  need  of 
the  standard  results.  Our  notation  follows  [  Parlett,  chap.8  ]. 

For  any  square  complex  matrix  A  and  an^  scalar  (  a  ]  ( called  the  shift )  that  excludes  A ’s  eigen¬ 
values,  the  associated  QR  transformation  A  A  is  defined  as  follows: 

i)  let  A  -  a  I  =  QR ,  the  unique  unitary  upper  triangular  decomposition  with  the  diagonal  ele¬ 
ments  of  R  being  positive. 

ii)  define  A  =RQ  +  al  =  Q*AQ . 

One  important  property  of  the  QR  transformation  is  that  both  the  upper  Hessenberg  form 
(A  =(fly)  with  aij  =  0  is  i>j+l )  and  the  Hermitian  form  ( A  =  (a^)  with  Oij  =a^- )  are  preserved.  Our 
concern  here  is  with  real  symmetric  tiidiagonal  matrices,  A  =  T,  and  this  form  is  preserved  in  the  QR 
transformation  since  T  is  both  upper  Hessenberg  and  Hermitian.  Only  real  shifts  are  considered  in  our 
investigations. 

2.2  Eigenvalue  deflation  and  eigenvector  calculation 

A  well  known  result  (  sec  [  Wilkinson,  pp  469-471  ]  )  connects  QR  with  eigenvalue  deflation  and 
eigenvector  computation. 

Definition  2.1:  A  symmetric  tridiagonal  matrix  T  is  called  unreduced  if  its  subdiagonal  elements  are 
nonzero. 

Remark:  Wlten  T  is  unreduced  the  QR  tran^ormcuion  is  well  defined  for  all  shifts  a  because  the  first 
n-1  columns  ofT  -  ol  are  linear  independent  for  all  o. 


Lemma  2.1:  ( QR  and  deflation ) 

Let  T  be  unreduced  and  f  be  the  QR  tran^orm  ofT  with  shift  a,  i.e. 
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f:  =  Q*TQ  =  RQ +al  (2.2.1) 

where  T  -  al  =  QR.If  a  =  X,  an  eigenvalue  of  T,  then 

(1)  last  row  off  has  the  form  (0 . 0,  X); 

(2)  last  column  of  Q ,  namely  q. ,  satisfies 

Tq,  =  q,X.  (2.2.2) 

Here  Q  is  orthogonal  and  R  is  upper  triangular. 

Since  f  =T  ©X  where  f  has  order  one  less  than  T,  we  say  that  X  has  been  deflated  from  r  in  one 
sweep  of  QR  transformation.  It  is  clear  that  the  spectrum  of  T  consists  of  the  remaining  eigenvalues  of  T. 
Also  from  Lemma  2.1,  we  see  that  when  X  is  deflated  fiom  T,  its  corresponding  eigenvecux*  is  revealed  in 
Q ,  namely  its  last  column  q„ . 

23  Some  examples 

In  this  subsection,  we  show,  by  example,  that  Lemma  2.1  is  not  a  reliable  guide  to  results  in  finite 
precision  computation.  Example  2.1  will  show  a  successful  deflation  and  Example  2.2  will  show  a  failure. 
Example  2.3  will  exhibit  the  success  of  deflation  cm  the  failed  case  in  Example  2.2  after  two  sweeps  of 
QR  have  been  applied.  Example  2.4  is  an  interesting  case  of  success  despite  having  a  shift  a  that  is  an 
exact  eigenvalue  of  several  of  T’s  leading  principal  submatrices. 

The  data  given  in  the  following  matrices  have  been  multiplied  by  10*  for  the  purpose  of  better 
presentation.  The  transformed  f  here  is  generated  by  the  numerical  implementation  of  QR  called  TQR, 
which  will  be  described  in  section  3  J!. 

Example  2.1:  ( the  successful  case ) 

6683.3333  14899.672 
14899.672  33336.632  34.640987 

34.640987  20.028014  11.832164 

11.832164  20.001858  10.141851 

10.141851  20.002287  7.5592896 
7.5592896  20.002859 

The  eigenvalues  of  this  matrix  are: 

Xi  =  0,  X2=10,  X3  =  20.  X4  =  30,  X5  =  40,  X«  =  40000. 

The  shift  is  Xi  =  0.  The  matrix  f  after  one  QR  sweep  is: 

39999.925  54.726511 
54.726511  33.404823  8.3017268 

8.3017268  24.730751  8.8065994 

8.8065994  21.646903  7.2175779 

7.2175779  20.292461  -7.943d-12 
-7.943d-12  -2.344d-15 

The  last  row  of  Ts  is  negligible  as  we  expected.  For  comparision,  here  is  computed  by  a  method 
other  than  TQR. 
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7*6 


•  39999.925  54.726511 
54.726511  33.404823  8.3017268 

8.3017268  24.730751  8.8065994 

8.8065994  21.646903  7.2175779 

12115119  20.292461  -1.113d-14 
-1.113d-14  9.52(W-13 


The  matrix  elements  of  these  two  transformed  are  almost  identical  except  the  bottom  ones.  How¬ 
ever,  they  are  negligible. 


Example  2.2:  f  The  failed  case  1 

The  matrix  T  is  the  same  as  the  one  ui  Example  2.1.  The  shift  is  X«.  The  matrix  f  after  one  QR 
sweep  is: 


19.989995  14.142133 
14.142133  20.003002  11.832160 

11.832160  20.001858  10.141851 

10.141851  20.002287  7.5593584 

7.5593584  20.730517  -17056153 
-17056153  39999.272 

J 


The  last  subdiagonal  element  is  not  negligible.  For  comparison,  here  is  fs  computed  by  a  method 
other  than  TQR. 


'  19.989995  14.142133 
14.142133  20.003002  11.832160 

11.832160  20.001858  10.141851 

10.141851  20.002287  7.5592896 

7.5592896  20.002859  -1.608d-13 
-1.608d-13  40000.000 

V,  J 


Examples  2.1,  2.2  have  shown  that  the  transformation  with  a  =  X,]  is  stable  and  the  one  with  a  = 
is  unstable.  Examination  of  the  eigenvalues  of  the  leading  principal  submatrices  of  7$  (  see  table  23.1  in 
Appendix  A )  reveals  that  matched  the  biggest  eigenvalues  of  7-^,1^,  and  T^lo  almost  full  working  pre¬ 
cision.  On  the  other  hand  is  not  close  to  any  eigenvalues  of  7-^,7 ^  and  7  5. 


Example  25: 

The  tiidiagonal  matrix  r«  is  the  same  as  that  in  Example  2.2,  We  applied  the  QR  transformation 
once  more  to  the "  7^'^  "  exhibited  in  Example  2.2  keeping  the  same  shift  =  40000.  The  resulting 
matrix  is: 


19.979990  14.142125 
14.142125  20.006003  11.832161 

11.832161  20.003716  10.141851 

10.141851  20.004574  7.5592897 

7.5592897  20.005717  8.425d-15 
8.425d-15  40000.000 
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The  last  subdiagonal  element  is  now  negligible. 

Example  2J  has  shown  that,  in  one  case  at  least,  TQR  will  take  two  sweeps  to  get  the  deflated  f 
Examination  of  the  eigenvalues  of  the  leading  principal  submatrices  of  obtained  by  first  QR  sweep 
with  a  =  reveals  that  the  first  QR  sweep  does  no  more  than  destroy  the  closeness  of  the  eigenvalues  of 
r«’s  leading  principal  submatrices,  ( see  table  23.2  in  Appendix  A ). 

Example  2.4;  ( successful  deflation  in  an  interesting  case ) 

The  matrix  T  in  this  example  is  the  well  known  second  difference  matrix.  The  data  are  the  original 
ones. 

'  -2  1 
1  -2  1 

Ts  =  1-2  1 

1  -2  1 
1  -2 

The  eigenvalues  of  this  matrix  are: 

Xi  =  -2-V3,  X2  =  -3,  h  =  -2,  X,  =  -l,  X5  =  -2  +  >/3. 

The  shift  is  X3  =  -  2.  The  matrix  fs  after  one  QR  sweep  is: 

‘  -2.0000000  1.4142136 

1.4142136  -2.0000000  0.70710678 
fj  =  0.70710678  -2.0000000  1.2247449 

1.2247449  -2.0000000  0.000000 

0.0000000  -2.00000000 

One  can  verify  that  X3  is  also  an  eigenvalue  of  the  first  and  the  third  leading  principal  submatrices  of 
7*5.  This  example  shows  that  even  if  the  shift  is  an  eigenvalue  of  some  of  the  leading  principal  sub¬ 
matrices,  nevertheless  QR  deflates  T  in  one  sweep. 

2.4  Spectral  properties  of  a  symmetric  tridiagonai  matrix  T 

We  give  here  several  results  that  we  need  later.  They  are  applications  of  Cauchy’s  Interlace 
Theorem,  sec  [  Cauchy,  vol.  2  ]. 

Theorem  2.1:  IfT  is  unreduced  then 

(1)  the  eigenvalues  of  T  are  distinct; 

(2)  the  eigenvalues  ofT's  consecutive  leading  principal  submatrices  interlace  with  each  other. 
For  proof,  see  [  Parlctt,  sec.  7-7,  sec.  7-10  ],  ( Cao,  chap.  2  ]. 

Definition  2.2:  spread ( T  )  =  ( T  )  -  (T). 

Theorem  2.2:  IfT  is  unreduced,  then  the  subdiagonal  element  P*  satisfies  the  following  inequality: 


if  n  >2,  I  P*  I  <  spread{T)l2, 
if  n  =  2,  I  P*  I  ^  spread  (  T  )  /  2. 


for  k  =1, ...,  n ; 


Definition  U:  Let 
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X*(<T):  =  det(rt  -  a/i ): 

5*;  =  X*(o)/X*-i(o)- 

where  is  the  k~th  leading  principal  subnuurix  ofT. 

Since  Xt  is  not  monic  it  differs  from  the  characteristic  polynomial  of  by  a  factor  ( -  1  )* . 
Notation  2.1:  V/e  denote  the  eigenvalues  ofT  by  A,-  and  label  them  so  that 

Aj  <  ^  ^  An . 

Notation  2.2:  Here  we  denote  the  bottom  element  of  the  normalized  eigenvector  ofXi  by  co, . 

In  our  applications  we  need  1  /  co,-  so  we  present  our  results  in  that  form. 

Lemma  2.2:  IfTis  unreduced,  then 


(2.4.1) 

(2.4.2) 


J_ 

(of 


xl(^) 

X-i(^)‘ 


(2.4.3) 


For  proof,  see  [  Parlett,  chap.7  ]. 


Theorem  2J:  Suppose  T  is  unreduced.  Let  us  denote  the  eigenvalues  ofT^-x  by  so  that 
f  <  |42  <  •  •  •  <  |i«-l 


then 


and 


Proof: 


Since 


>  i 


(ot 


<  i 


Wj 


A2~  Ai 

Vii-h’ 

Aj  ~  Aj-i  A,>}  —  Aj 

Ai-M.-,  H.-A, 

K  •  ^n-l 

An  -  Ai 

\lx-h’ 

Aj  —  Aj  An  —  Aj 


Wi 


Aj  -Ui_i  iij  -Aj  ’ 
An-Ai 

An  -iU-1  ’ 

xU>h) 

X— 1(^  )’ 


i  =  1; 
i  *  I,  n; 
i  =«; 

1=1; 

/■  *  l,n; 
i  =  n. 

by  (2.4.3). 

II -1 


Xnio)  =  n(^y-®).  x-.-i(<y)  =  n(M^,-<T). 

/»!  /=! 


(2.4.4) 


(2.4.5) 


then 
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And 


/■i.y  *i 


n-l 


n  (^/-^) 

/«i.>  ■»« _ 

n(ny-:^) 


By  Theorem  2.1 


/-I 


1  =  1; 


i-I 


n 

y-i 


^2 _ h-  n  ^ 

^  /-^i  V-j-K 


i  *  l,n; 


n? — -> 


i  =n. 


Xj  <  Hj  <  Xj^.t,  j  =  l . n-1. 

Therefore  each  factor  in  the  products  in  (2.4.7)  is  positive  and  is  bigger  than  one. 
formulae  in  (2.4.7)  satisfy: 


lii-Xi 

X|  ~  X,-|  X,>|  —  X,- 
X,-|i,-i  M,-X, 

x»  -  X,-i 

X,-H.-,' 


i  =  1; 
i  *  l,n; 
i  =  n. 


The  rervganization  of  (2.4.7)  reveals  that 


n  htizh.  =  hizh.  n  h—h 

yVi  liy  -  Xi  Hi  -  X,  ]lj  -  X,  ’ 

\  -_Xj  yj  ^ 
yVi  X,  -nj  \ij  -Xi 


=  ( n  hzltL ) huh  (  fr  hj-\ 

x,  -p.-i  Xi-]ij  ’  \ii-Xi  jiL\ij-Xi 


^  K  _Xy  ^  Tj  b 

K-V^-l  /.I  X, -Py 


(2.4.6) 


(2.4.7) 


That  is  to  say  the 


i  *  l,n. 


Since  Xy  <  liy  <  Xy^.],  each  factor  in  the  above  products  is  positive  and  each  factor  in  the  pro¬ 
ducts  on  the  right  hand  sides  of  equal  signs  is  smaller  than  one.  That  is  to  say  the  above  formulae 
satisfy: 
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<  i 


~  ^  \ 


M,  -^ 


^  ^1 
Xn  -M»-l 


i  =  1 
i  ^  l,n. 
i  =  n 


3.  Implementation  of  the  QR  transformation 


This  section  develops  the  usual  implementation  of  the  QR  algorithm  applied  to  a  symmetric  tridiago¬ 
nal  matrix  7.  In  the  process,  the  important  intermediate  quantities  are  introduced  along  with  the  relations 
among  them. 

Most  of  the  material  is  standard,  see  [  Wilkinson,  chap.8  ],  [  Stewart,  chap.7  ],  [  Parlett,  chap.8  ]  and 
[  Golub  &  Van  Loan,  chap.7  ].  However  all  the  results  are  needed  in  the  next  section. 

In  the  following  discussion,  we  assume  that  all  the  tridiagonal  matrices  in  question  are  unreduced 
since  otherwise  the  problem  decouples.  We  also  assume  that  the  offdiagonal  elements  <n  )oiT 

are  positive. 


3.1  QR  factorization  of  7  -a/ 

The  desired  QR  decomposition  can  be  carried  out  by  pre-multiplying  the  tridiagonal  matrix  T  -al 
by  a  sequence  of  plane  rotation  matrices  R^  £n  )  defined  as  follows. 


r  1 


/ 


1 


Ck 

-Sk  C* 


k~lh  row 


(3.1.1) 


The  duty  of  ( 2  ^  k  ^  n  )  is  to  annihilate  the  (  k,k-l )  position  of  the  matrix  on  the  way  to  an 
upper  triangular  form.  The  formulae  in  step  (k)  are  important  for  the  analysis  in  section  4. 

Let 


7  -  <T/ 


tti-O  P2 

P2  “2-0  ■ 

■  • 

P«  cu-o 


It  can  be  shown  that  at  step  (k)  ( k  <  n  ); 
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where 


^k^k-l  ' '  RtiT  -  csl )  - 


C*  X 
Ck  Pi+l 

Pt+l  ■ 


%k  =  ( Jt*-1  +  P*  Ci  =  Ck-iCk^k  +  Jt(  Ct*  -  a ) 

Ck=nk-Ak  5*=Pik/|t 

Ji*  =  -J*p*ct_i  +  ct(  a*  -  a  ). 


At  step  (n): 


with 


I  ^2  C2  X 


R,R,.i-RjiT-an 


^2  C2  X 


(T-al)  =  Q 


Q  =  R'jR'i-Rl. 


We  collect  some  direct  consequences  of  this  decomposition  that  are  used  later. 


a): 

Xi(a)  =  a,  -  a  =  71], 

X*(cr)  =  2<k^. 

(II):  IfT  is  unreduced,  then 

0.  St  *  0.  2^k^n. 

(in):  IfT  is  unreduced,  then 

7tt(a)  =  0  ifandort/y  if  X*(®)  =  0,  1  <it  ^n. 

In  words  vanishes  only  at  eigenvalues  ofTt .  So  does  Cn+j  by  its  definition. 


(3.1.2) 


(3.1.3) 


(3.1.4) 


(3.1.5) 

(3.1.6) 

(3.1.7) 
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Lemma  3.1:  The  detailed  structure  of  matrix  Q  is: 


C,C2  Ci(-S2)C3  Ci(-J2X-53)C4 
$2  ^2^3  ^  3)^  4 

Si  CiC^ 

S* 


Ci(-52)-(-5«)  1 


it'll  C*— 1(  •^») 

■T.  Cm 


(3.1.8) 


3.2  Thei2^  transfwmation 

In  this  subsection,  we  look  at  the  inner  loop  of  the  QX  transformation.  The  formulae 
t  =  Q‘TQ  =RQ  +<sl  in  (2.2.1),  R  in  (3.1.4)  and  Q  =/?2  •  •  ■  in  (3.1.8)  suggest  a  way  to  transform  T 

intoT  without  forming  |2  • 

I^t  let  us  collect  the  essential  quantities  in  the  QR  factorization  derived  above.  If  we  define 

Ji  =  0,  Cl  =1,  Jii  =  tti  -<j 

then  the  necessary  elements  of  R^.,  l^k^n,  canbe  generated  by  the  following  loop: 

For  k  =  2, ... ,  n 

^*=(Jt*"-i  +  p*^)‘'" 

Ck  ~  ^k-l  /  %k 
Sk-^k  I  %k 

-  «*  =  -  «*P*C*_i  +  C*(  O*  -O). 

f 

Now,  let  us  look  at  the  diagonal  and  subdiagonal  elements  of  matrix  RQ  +al  to  find  out  the  formu¬ 
lae  needed  to  compute  the  matrix  elements  of  f. 

4 

By  the  formulae  (3.1.4),  (3.1.8),  RQ  +  al  can  be  presented  as  follows: 


»  <4 

'  C2  -S2Ci  •  •  •  Ci(-S2)...(-S,)  ■ 

52  C2C3  •  •  • 

53  •  • 

%m  C, 

Cm—ICm  Cp,_i(~5, ) 

b  4 

Sm  C, 

1.  4 

If  we  denote  the  diagonal  and  subdiagonal  elements  of  f  by  and  respectively,  direct  calculations 
reveal  that,  for  2  £  k  n , 

P*-l  =  %k^k-li 

d*_i  =  4*c*_ic*  +  +  o. 

=  ^c*_iC*  +c*_,c*P*st  +St*(a*  -a)  +  o,  C*  =c*_ic*p*  +54(0*  -o) by  (3.1.3), 

=  c*_iJtt_i  +  c*.iC*Pt j*  - c*^  a*  - o)  +  a*,  n*.,  =  ^*Ck  by  (3.1.3), 

=  ct_iit*_i-c*(ct(at  -a)-p*j*c*_i)  +  cq, 

=  “  ctTt*  +  a*. 


=  ^mSmi 

a,  =  n,c,  +0. 
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To  organize  the  computation  it  is  convenient  to  introduce  a  new  quantity  yt : 

Definition  3.1: 

The  relation  between  Yt  and  y^^.j  in  toms  of  and  is 

Yi  =  "iCi 

Y*  =  c^(ctt-a)  -  ytW.  2^*  Sn. 

The  above  fcmnula  is  derived  as  follows 

Jt*c*  =  ( -  + 1*(  a*  -  c  ))  c*  by  definition  of  jt*  in  (3.1.3), 

=  -^*P*ctc*_i  +  c/(a*-c), 

=  -  Jt  +  ct^  Ot  -  o  ),  since  Ct  Pt  =  "i-iJt  from  (3.1.3). 

Linking  all  the  above  relations  together,  a  compact  algorithm  emerges.  We  list  this  algorithm  in 
detail  in  (4.1.1.)  and  name  it  in  this  study  by  TQR.  It  is  essentially  the  algorithm  used  in  the  well  known 
EISPACK  collection  of  routines. 

Example:  In  the  following,  the  computed  ct’s,  tCt’s  from  TQR  in  Example  2  J  (  see  section  2  J  )  are 
listed.  The  resulting  T  is  exhibited  in  (2.32).  For  comparison,  the  correct  ct 's  and  nt  's  are  also 
listed.  The  resulting  f  is  exhibited  in  (2.3.3). 


computed  correct 


k  c*  Ct 


1 

+1.00000(x)000000d400 

+1.0000000000000d-t00 

2 

-9.1287087206375d-0l 

-9.1287087206375d-01 

3 

+7.9096456588243d-O4 

+7.90964565943(XW-O4 

4 

-2.338830021282W-07 

-2.3408762727625d-07 

5 

-8.0658942646820d-07 

+5.9381746504385d-ll 

6 

+4.2662106420853d-03 

-1.1223688023196d-14 

k 

n* 

1 

-3.3316666666667d+00 

-3.3316666666667d+00 

2 

+2.7399802074030d-06 

+2.7399802074446<i-06 

3 

-2.76734 19903274d-10 

-2.7697632729029d-10 

4 

-8.1803099953432d-10 

+6.0232682537306d-14 

5 

+3.2249815432264d-06 

-1.9058339689554d-17 

6 

-1.705630831781  W-02 

-1.6080662764449d-17 

3.3  Analysis  of  the  QR  transformation 

Section  3.2  reveals  the  relations  between  the  quantities  of  Ct,St,  and  the  matrix  elements  ctt ,  pt 
in  one  step.  This  is  inadequate  if  the  analysis  in  terms  of  a  is  needed.  In  the  next  lemma,  we  present  several 
matrix-vector  relations  b^ween  all  the  intermediate  quantities  generated  in  the  QR  process  and  the  matrix 
T.  It  is  these  relations  which  will  help  us  understand  QR  more  deeply.  These  relations  also  tell  us  the 
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stracture  of  an  eigenvector  when  a  happens  to  be  an  eigenvalue  of  T. 

Recall  the  partial  reduction  ofT-al  to  upper  triangular  form  as  it  appears  at  step  (k).  It  is  given  in 

ai^). 

The  product  of  the  plane  rotation  matrices  R2>  —  >  satisfies: 

C2  -S2C3  •  *  '  ' 

S2  C2C3  •  •  • 

S3  •  • 


Ryti 


Ck-lCk  C*_i(-St) 

Sk  Ck 


*11— k 


(3.3.1) 


Thek-th  column  of  (3.3.1)  plays  a  crucial  role  in  our  analysis. 

Definition  3.2:  We  denote  by  jk  the  following  vector  in  R* . 

Ci(-S2)  •  •  •  (-St)  1 


y*:  = 


C2(-S3)  •  •  •  (Sk) 

Ck-l(-Sk) 

Ck 


(3.3.2) 


Definition  3.3:  We  denote  by  y*  the  following  vector  in  R" . 

Sk-  -  ^2  ■  •  ‘^t  «t  = 

Equating  thejfc-./:  row  of  vJ-1-2),  the  following  matrix- vector  relations  emerge. 


y* 

0 

0 


(3.3.3) 


Lemma  3.2:  If  T  is  a  symmetric  tridiagonal  matrix  of  order  n  and  a  is  an  arbitrary  real  number,  then 
the  quantities  derived  up  to  step  k  in  TQR  sati^ 


Tik  -  h<y  =  «tet  +  CiPtt,et+,. 

k  <n. 

(3.3.4) 

fl  Tfk  -  yt<y  II*  =  ni  +  c*^pt*+i. 

k  <n. 

(3.3.5) 

(SkYiTyt-aSk  )  =  JitCt  =  Yt  =  (SkYTh 

-  0,  A  S  n , 

(3.3.6) 

Tfn  -  y»<y  = 

(3.3.7) 

Proof:  Equate  the  k-th  row  on  each  side  of  (3.1.2)  and  transpose  to  get 

(r-a/ K^2  •  ■  •^t  )et  =  (T-al)yk  =  rtt** +i^iPt+iet+i. 

Since  (3.1.2)  holds  fork  <n,  (3.3.4)  is  true  for  ifc  <  n .  (3.3.5),  (3.3.6)  are  the  direct  results  of  (3.3.4). 
(3.3.7)  is  a  special  case  of  (3.3.4)  since  P,+i  =  0  when  k  =n.  ■ 

We  use  notation  y^  instead  of  qt  since  the  former  is  slightly  different  from  the  k  -th  column  of 
matrix  Q .  When  i  =  n ,  y,  =  q,. 
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Lemma  3 J  used  T  and  f], .  Equally  important  is  the  relation  between  Ti^ ,  the  leading  kxk  submatrix 
of  r.andyt. 


Coronary  1  of  Lemma  32:  For  any  real  a, 

TkJk  -  y*<J  =  l^k^n.  (3.3.8) 

(y*)'(7*y*  -  ay*  )  =  Jt*ct  =  Y*  =  (y*)'7'ty*  -  (3.3.9) 

Proof:  Taking  the  first  k  rows  of  (3.3.4)  and  using  the  notation  in  (3.3.2),  the  formula  in  (3.3.8)  is 
obtained  fork  <n.  When  k  =  n  it  is  the  case  in  (3.3.7). 

(3.3.9)  is  the  direct  results  of  (3.3.8)  and  (3.32).  ■ 


Coronary  2  of  Lemma  32:  If  a  is  an  eigenvalue  ofTk .  then 

Tkyk  -  y*cf  =  0. 

TSk  -  yt<T  =  c*P*+ie*.^i.  II  Ty*  -  y*c  II  =  I  c*p*+i  I, 
a  =  (y*yT*y*  =  (y*)‘ry*. 

Direct  consequences  of  (3.1.7),  (3.3.8),  (3.3.4),  (3.3.5),  (3.3.6)  and  (3.3.9). 

Definition  3.4:  We  denote  by  t*  the  following  vector  in 

_  r  n*  ■ 

»  a 


(3.3.10) 

(3.3.11) 

(3.3.12) 


(3.3.13) 


4.  Properties  of  Cj^  and 

The  main  results  are  (4.1.12)  and  the  tight  bounds  on  the  doivatives  of  n*,c*  and  llt*ll  with  respect 
to  a.  See  Theorems  4.1-4.5. 

4.1  Basic  properties. 

We  label  the  eigenvalues  of  T*  so  that 

<  •  •  •  < 

For  easy  reference,  TQR  and  relations  (3.1.7),  (3.3.8),  (3.3.2)  are  restated  here. 


TQR:  (4.1.1) 

Ji  =  0,  Cl  =  1,  iti  =  a, -a,  Yi  *  TiiCi  =  Wi 

for  k  =  2, ... ,  n  do 

=  P*  / 

C*  =  lt*-l  /  kk 

n*  =-5*PikC*_i  +  c*(ot*  -  a) 

Yjk  =c*^a*-a)-j*^*_, 

^-i=Y*-i~Y*  +  at 
~  Pt-I  =  kk^k-l 

a,=Y«+a. 
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Relation  (3.1.7):  If  T  is  unreduced,  then 

=  0  ifandort/y  if  <J  = 
c*(c)  =  0  ifandon/>  if  a  = 

Relation (3 J.8):  (I*  -a/t)y*  =  it*  e^*^ 

Relation  (3  J.2): 

Ci(-S^  •  •  •  (Sk) 

y*  = 

Ck-l(-Sk) 

Ck 


Some  preliminary  results: 


I  ^i  ^k, 

l^i^k-l,  2^k^n. 
l^k  ^n. 


yLyk  =  1- 


Corollary  of  Relation  (3  J.8):  When  <r  =  Xj^*\  then  y*  i«  (4.1.4)  is  its  eigenvector  for  Tk . 


(4.1.2) 

(4.1.3) 

(4.1.4) 


(4.1.5) 


Notation  4.1;  We  denote  the  bottom  element  of  the  normalized  eigenvector  of  by  co.  jt, 
i  =  1.2, ....  A,  ik  =  1.  ....n. 


Corollary  of  Relation  (3  J.2):  When  a  - 

(Oijt  =  ct  =  CkiU^^)  *  0  i  =  l,2 . ik,  i:  =  l,2 . n.  (4.1.6) 

This  is  the  direct  consequence  of  (4.1.S)  and  (4.1.3). 


We  now  discuss  the  smoothness  of  «* (o),  c*(a),  for  o  e  ( - »®,  <» ). 


Lemma  4.1;  IfT  is  unreduced  and  Ck.Sk.Hk  are  the  functions  computed  by  TQR  (  see  (4.1.1) ),  then 
Ckt  SktKkt  llfk  II  are  real  analytic  on  R. 

Proof:  It  may  be  verified  that 

re*  =  det*[  r*  -  a/fc  ]/det[ (T*.!  -  olk-i)^  +  Pkek-i*k-i]> 

Therefore  ni  is  a  rational  function  of  order  2kl^-2)  with  no  poles  on  the  real  axis.  Thoefore  tc^ 
and  4k  =  (J^k-i  +  Pk  )*^  are  analytic  on  R  for  Jk  =  2, ....  n .  Since  4k  >  0.  H  tk  H  >  0  it  follows  that 
Ct.Jt  and  lit*  II  are  also  analytic  on  R.  Note  that  n*-»o^  as  a-*oo. 

We  want  to  show  that  n*  is  a  weighted  harmonic  mean  of  the  ( -  a  )^,  i  =  1. ....  k . 


Premultiply  both  sides  of  (4.1.4)  by  (  T*  -  <r/*  )~*Jtk  to  find 

-^  =  ( r*  -  a/k  )-*«^*\  (or  a  ^ 

«k 

A  consequence  of  the  spectral  factorization  is 

(e*»>)'(r*-a/*)-^ei*>  =  for 

iMiin)  ’-ay 

Since  y*y*  =  1  by  (4.1.5),  (4.1.7)  yields 

-T  =  (ei*>y(r*-a/*)-"ei*>. 


"k 

Combine  (4.1.8)  and  (4.1.9)  to  find 


(4.1.7) 


(4.1.8) 


(4.1.9) 


1 


W.4 


nl  ■  W*>-<T 


for  o  * 


(4.1.10)B 
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Applyin^  Lemma  4.1,  a  fundamental  relation  between  Ck  and  their  derivatives  with  respect  to  a 
is  obtain^.  /  (a)  =  df  (a)/da. 

Lemma  4.2:  For  all  real  a 

lik  Ck  -  c*  =  1.  (4.1.11) 

Proof:  Differentiate  (4.1.4) 

-7k  +  (Tk-<ifk)yk  = 

Multiply  by  ( y^  )'  and  recall  the  definition  of  yk  in  (4.1.S) 

-(y*  )'y*  +  (y*  )'(7t-CT/»  )yi  =  ctn*. 

The  result  follows  since  (  y*  )'y*  =  1  and  (  y*  )'(  T*  -  a/*  )  =  «*(  ei*^ )'  by  (4.1.4).  ■ 


Recall  that  t*:  =  (Jtt,c*Pt+i)‘. 

Corollary  1  of  Lemma  4  J:  For  any  real  c 

fittll  lift  11  tsin/:(tt.U)l  (4.1.12) 

=  I  t*x  tt  I 

=  I  >t*c*P*+i  “  ^*CtPt+i  1 
=  I  I  by  Lemma  4.2 

=  Pt^f  ■ 

Discussion; 

Relation  (4.1.12)  tells  us  that  fitkil  wiU  be  huge  if  11(^11  I  sin  (tk,tk)l  is  tiny  and  P^.,.!  is 
moderate.  In  cnder  to  understand  the  behavior  of  lift  II  deq)ly,  a  detailed  analysis  has  be^  made  on  the 
behaviors  of  Kk  and  Ck ,  of  which  II  II  consists. 

4.2  Properties  of 

In  this  subsection  we  investigate  the  behavior  of  re* ,  2  ^  /fc  S  n  for  a  e  (  -  oo,  <» ). 

Corollary  2  of  Lemma  4.2: 

)  =  -  -TTwx  =  -  1  ^ ‘  ^ 

c*(VO  “'•* 

Proof:  Sincen*(X/*’)  =  Oby (4.1.2), relation (4.1.11) at  becomesct(  V*^)jtt(X.i^*^)  =  -1. Since 
c*  ( )  =  Wi  *  ^  0  by  (4. 1 .6),  (4.2. 1)  is  obtained.  ■ 

Corollary  3  of  Lemma  42:  For  any  eigenvalue  ofTk , 

itk(o)  =  Jti(X,<*>)(o-X<^*))yi(a)y*(X,«*^). 

Proof:  Rewrite  (4.1 .4)  as 

(  r*-a  +  a-X<<*>  )y*(X/‘>)  =  e*nt(Xi(*))  =  0. 

Premultiply  by  yt(  o )  to  get 

n*(  a )  ct(  Xi<‘> )  =  -  «T  -  X/*) )  yi(a)y*(Xi(*)). 
and  then  apply  Corollary  2  of  Lemma  4.2. 
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CoroUary  4  of  Lemma  A2'. 

=  0.  (4.2.2) 

Proof:  Differentiate  (4.1.11):  jt*  c*”  -  c*  it*  =  0.  Set  o  =  X/*\  then  jt*  =  0  by  (4.1.2)  and  c*  =  co/jt  *  0 
by  (4.1.6).  ■ 


We  now  turn  to  the  behaviour  of  tc^  for  all  other  values  of  a. 
Lemma  43t 

nkiti  <  0,  for 
Proof:  Differentiate  both  sides  of  (4.1.10) 


for  o  ^ 


Differentiate  right  hand  side  of  (4.2.4) 


3  y - - 


Differentiate  left  hand  side  of  (4.2.4) 


Therefore  fw  c  ^ 

1  «  3  i, 

*s  “  -a  " 

*  ».  .  . 

k 

=  3  1 


if* 


1  ^  3  /  ^  v2 

--T^/k  +  —(«*)• 
tt* 


3  /,  r  \2 


(4.2.3) 


(4.2.4) 


(4.2.5) 


(4.2.6) 


«*  K*  .t'i( 

w?.* 


(jt*r 


T  _  3  r  _L  12 

.•t'.(X.^*^-a)^■t'.(X.^‘)-a)^ 


by  (4.2.5)  and  (4.2.6) 

using  (4.1.10) 


=  3[  K 


‘  “a  ^2  i  “ijt 


-(L 


* 


S  0. 


by  Cauchy-Schwarz  Inequality. 


)^  ]  using  (4.2.4) 


Moreover  equality  holds  if  and  only  if  the  following  two  vectors 


zi  =  ( 


©u 


©t.* 


)' 


Z2  =  ( 

are  proportional.  That  is  to  say: 


©u 


©M 


or 


(X|*)-<T)^’  ’(Xi*>-a)^ 

-  a  =  •  •  •  =  X^*^  -  a, 

X.<*’  =  •  ■  ■  =  xi*). 


)' 


That  contradicts  the  conclusion  of  Theorem  2.1  in  section  2.4.  Therefore 
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n* 

Recall  from  (4. 1 .2)  that  «*  (a)  ?6  0  for  <r  ^  X/*\  Therefore 

-  JT*  jt*  >  0.  for  o  ^  ■ 

Corollary  of  Lemma  4  J: 

Kt  *  0.  for  (4.2.7) 

Lemma  4J  shows  that  the  algebraic  function  nt(<T)  is  like  the  characteristic  polynomial  of  in  that 
it  vanishes  at  the  eigenvalues  of  T^.  Moreover  it  is  altematingly  concave  upwaid  and  downward  in  the 
intervals  bounded  by  the  eigenvalues  of  T/^. 

The  next  result  is  a  direct  corollary  of  the  proceeding  lemmas  and  Theorem  23  in  section  2.4.  We 
illustrate  it  in  Fig.  4.1  which  shows  that  attains  its  extreme  values  at  the  i  -  1 . k . 


Theorem  4.1:  The  derivative  of  the  function  Jt*  (a)  computed  by  TQR  satisfies 


[  «*( 


(4.2.8a) 

I  l.k,  (4.2.8b) 

(4.2.8c) 


and 

I  «;(!.,<*>)  I,  ae  (-00,  A./*)]; 

I  n'tia)  I  5  max  [  I  )  I,  I  iii(  )  N  ®  ^  ],  i  =  1 . k-l; 

ln;(Xi*>)l.  ae[Xi*>.  oo). 
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Proof:  By  (4.2.7):  n*  ^  0  for  a*  Apply  Corollary  4  of  Lemma  4.2 

jc*  =  0  for  a  =  V*^- 

Thus  on  each  [  X/JJ  ]  for  i  =  1 . k-l,  w*  has  its  only  stationary  points  at  the  ends.  For  the 

intCTvals  (  -  o®,  Xf*^  ]  and  [  Xi*\  <» )  the  stationary  points  are  X/*^  and  Xi*\  Therefore 

I  ni(o)  I  ^  max  [  1  TC|[(X/*>)  I,  I  Jti(X/t{ )  I  ],  for  t  =  1,  ...,fc-l. 

Use  (4.2.1)  to  find  that 

tt*(  ^**’ )  =  -  1  /  tOjjt.  for  i  =  1 . k , 

Apply  the  bounds  for  1  /  oi^ji  in  Theorem  23  on  Ti^,  the  formulae  in  (4.2.8a),  (4.2.8b)  and  (4.2.8c) 
are  revealed.  ■ 


We  now  give  a  pointwise  bound  on  that  reveals  the  role  of  the  distance  of  a  from  the  spectrum  of 
Tt  and  firom  the  spectrum  of  Tt-i.  Some  preliminary  results  are  restated  for  easy  reference. 


(A) 

t  =  (7t*2-l  +  P*^)*'" 

see  definition  in  (4.1.1) 

(B) 

see  definition  in  (4.1.1) 

(C) 

see  definition  in  (4.1.1) 

(D) 

,  St  / 

c*  =  ^"*-1 

derivative  of  (B) 

(E) 

’tjC*  —  C*7t*  =  1 

for  any  real  a 

(from  (4.1.11)) 

(.F) 

(r*  -  <J/*  )y*  =  n*  e^*> 

for  any  real  a 

(  from  (4.1.4) ) 

Definition 4.1:  p.*:  =  M*(<*)  =  min  IX/*)-<Jl,  X  =  l,2,...,n. 

lii  ik 


\i*  <  n* 


(G) 

Proof  of  (G):  (F )'  (F )  yields 


for  CT  *  X/*> 

(y*  )'  (T*  -  a/*  )^yt  =  nl 


Then 


=  min  (o-X,^>)2 

ISi  ik 

=  Knniin-ahf) 

<  (yknn  -  a/*)2y* 


by  Definition  4.1, 

y*  cannot  be  an  eigenvector  since  ct  *  X/*^ 


We  now  present  one  of  our  key  technical  lemmas. 

Lemma  4.4:  For  any  real  a, 

Pt(a)  I  I  =  I  Jt*  I  CT  =  X,^*) 

p*(o)  I  jt*  I  <  1 71*  1  a  Xif*> 

H*_i(o)  1 7t*  I  <  s*^  I  It*  1  +  (pt^_i(a)  +  p**)*'^ 


(4.2.9a) 

(4.2.9b) 

(4.2.10) 


Proof  of  (4.2.9a, b):  When  ct=  X/*^  then  tc*  =  0  by  (4.1.2),  (i*  =  0  by  definition  and  (4.2.9a)  holds.  Now 
suppose  that  a  *  X/*\  Premultiply  (F)  by  ( 7*  -  a  /*  )"'7it"'. 


=  (T*-a/*)''ei*>, 

^k 


(4.2.11) 


Differentiate  (4.2.11), 
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Jk  >Ct  -  yt  rtt 


=  (r*  -  a/*r"ei‘\ 


(4.2.12) 


Differ^tiate  y*yt  =  1  to  obtain  y*y*  =  0.  ftiemultiply  (4.2.12)  by  y*,  getting 


=  (y*)‘(7’t-<T/*)-‘yt/7tt. 


since  yiyi  =  0 
by  (4.2.11). 


Thus 


and 


-«*  =  (y*)'(7*-<TAr‘y*nt. 


(4.2.13) 


l«;  I  =  l(y*)'(7’*-o/t)-'yt  I  l7t*  I. 

Since  a  ^  X/*\  y*  will  not  be  an  eigenvector.  Therefore 

I  (y* )'( n-<jh )'‘  y*  I  <  nm  I  v'  ( r*  -  c/* )-‘  v  i 

=  ll(r*-a/*)-‘ 

=  1  ^  J_ 

min  I  X/*'  -cl  14 

1  iiik 

Put  this  inequality  into  (4.2.13)  and  multiply  by  to  obtain  (4.2.9b).  ■ 

Proof  of  (4  J.  10):  When  o  =  then  and  (4.2.10)  is  immediate  since  p*  ^Oandj**  Itt*  I  ;*0 

by  (3.1.6)  and  (4.1,2).  , 

Now  suppose 

I  jti  c*  +  1  I  =  I  c*  re*  I 

=  I  «*_:  re*  I 

%k  lit-i 
Sk  I  C*rt*  I 


by  (E) 
by  (D) 

by  (4.2.9b) 

by(B) 


Hence 


.  ,  Sit  I  c*re*  I  , 

I  c*  I  <  - +  1. 

14-1 

Since  a  *■  V*“*\  c*  0  by  (4. 1 .3).  Then  the  above  inequality  can  be  rearranged  as 


I  rei  I  <  j*^  I  re*  I  + 
=  j*2  I  re*  I  +  p*_i 
=  j*2  I  re*  I  +  tt*_, 


14-1 
I  c*  I 

\k 

I  re*^.,  I 


I  It*-,  I 


by  (B) 
by  (A) 


.1/2 


=  s**  I  re*  I  +  (  \ii.i  +  ) 

<  s*^  I  re*  I  +  (  n*^_,  +  p*2  )>'*. 


by  (G) 
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The  geometric  interpretation  of  inequality  (4.2.9b)  is  illustrated  in  Fig.  4.1.  Since  the  graph  of  is 
concave  downward  or  concave  upward,  we  have 


,  -  ,  I  It*  I  I  I 
'  «* '  =  — r—  < 


for  a* 


The  die  is  defined  on  Fig.  4.1. 
Theorem  4.2:  For  any  real  a 


.  ispreadiTie)  3,,  spread  (Tt)  ^  , 


(4.2.14) 


Proof;  Recall  the  definition  of  TCjt  in  (4.1.1):  =  c*(a*-cy)  -  |3*JtCt_i.  Thus 

I  Jt*  I  =  I  c*  (ot*  -a)  -  P*s*c*_,  I. 

<  I  a*  -  a  I  +  I  pt  1 ,  since  I  c*  1  ^  1  for  2<k  <n 

=  I  <Xt  -  -  d  I  +  I  p*  I ,  choose  closest  eigenvalue 

S  I  I  +  (1*  +  spread  (Tit)  1 2,  by  Theorem  2.2 

<  spread(Tit)  +  14+  spread (Tt)  1 2,  since  a*  e 

=  4*  +  3  spread(  T*  )  /  2. 

Now  (4.2.9b)  in  Lemma  4.4  yields 


,  ,  ,  ^  I  Ttt  I  ,  3spread(Ti) 
I  Jt*  I  <1  -  <  1  + 


lit 


2m* 


(4.2.15) 


Next  we  use  the  definition  of  in  a  different  way. 

I  I  <  I  c*  I  I  a*  -  a  I  +  I  P*  I. 

=  I  cj  I  I  a*  -  -  a  I  +  I  Pt  I ,  choose  closest  eigenvalue; 

5  I  Ct  {(spread (Tit)  +  4t-i )  +  spread (T^) 1 2,  since  X)*"*^€  [X/*\X;J*M- 

So 

s^  I  TCt  I  <  st^  I  Ct  {(  spread (Tt)  +  4t-i  )  +  s^  spread (T^)  /  2 

S  (spread (Tit)  +  4t-i)/2  +  spread(Tt)  1 2,  since  I  St^t  >  ^  1/2. 

=  4t-i/2  +  spreadiJit).  (4.2.16) 

Next  we  bound  the  term  (p-t-i  +  P*  )'^  on  the  right  hand  side  of  (4.2.10). 


(m"-i  +  P*)"^  ^  +  (spread (Tit)  1 2)^  <  4t-i  +  spread (T,)  /  2.  (4.2.17) 


Add  (4.2.16)  and  (4.2.17)  to  get 

St"  I  Jtt  I  +  (4t^-,  +  Pt")*'"  <  f  (4t-i  +  spread  (Tit) ).  (4.2.18) 

Substitute  (4.2.18)  into  (4.2.10)  to  find 
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I.;,  <  ). 

2  M*-i 

The  combination  of  (42.15)  and  (4.2.19)  yields  (4.2.14).  ■ 


(4.2.19) 


Note  that  and  I4_i  cannot  vanish  simultaneously  by  Theorem  2.1  in  section  2.4. 

Theorem  4.2  shows  that  nt(cT)  can  only  be  huge  if  a  is  very  close  to  an  eigenvalue  of  and  an 
eigenvalue  of  r^-].  For  instance  when  I  1^1/'^,  (4.2.14)  implies 

M*  ^  3^6  14-1  ^  3^ 

spread(Tt)  ^  2(  l-'^)’  spread(Tt)  ^  2-3'^' 

43  Properties  of  ct^.]  and  St+t 

In  this  subsection  we  investigate  the  behavior  of  ct^-i,  1  ^  A:  ^  n-1  as  functions  of  the  shift  a. 
For  easy  reference,  we  collect  the  previous  results  that  are  needed. 


/tf\  '  4+1^*+!  ' 

4+1  ~  ~  6 

^+1^ 

,,,  »  3  Jt*  Pt+l  >  '  .2  _ _ 

(/)  ct+1  =  — -j —  iiitr  +  -3—  lit. 

S*+l  S*+l 


Pt+l 


(J)  14  K*  <  0  for 

(K)  «;(X,<*))  =  --f- 

Olijt 

Lemma  4  J:  Far  a  *  /  =  1 , ....  A: , 


Proof: 


4+1  *  0. 


-  '  _  3  Ttt  P*+l  /  \2  . 

4+1  ~  ,j  '^*^"^63 
s*+i  4+1 


derivative  of  (Q 

derivative  of  (D) 

( from  (4.2.3) ) 

(  from  (4.2.1)  ) 


by  (I). 


(4.3.1) 


However 


I  4+1  I  = 


3  I  It*  I  Pt+i  ^  ^2  .  Pt+i 


^*+1 


(«*)“'  + 


^+1 


I  71*  I, 


by(J) 


Lemma  4.6: 


> 

^'+1 

>  0. 


I  Tt*  I. 


4'+i  =  0. 


by  (J)  again. 


1  ^  ^  il. 


(4.3.2) 


Proof:  Use  (I)  and  observe  that  tc*”  ( )  =  0, 1  S  /  <  by  Corollary  4  of  Lemma  4.2  and 
7t*(  )  =  0, 1  S  i  5  A:  by  (4.1.2).  ■ 


Lemma  4.7:  c*4i  ( a )  teaches  its  extremes  at  the  and  then 

«t(  V*’ )  1 


4+1  (^‘*')  = 


Pt+i 


P*+i  ©ij 


1  S  i  5  A 


(4.3.3) 


Proof:  Since  c*Vi  €  C'(  -  oe,  00 ),  Lemma  4.5  and  Lemma  4.6  show  that  ci+i  teaches  its  extrema  at 
for  i  =  1 .....  A .  Therefore 

4+1  ( )  =  r2--^i-^.{,2  ’'*(  V*’ ).  by  (D)  and  (A) 

( ^*  +  P*+i ) 
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1 

Pt+i  Wijk 


by  (4.1.2)  and  (4.1.3) 
by  (K)  ■ 


Theorem  4.3:  The  derivative  of  function  ct+i(CT)  computed  by  TQR  satisfies 


<  [Pt+ictVi(Xt<*))]"  < 


(4.3.4a) 


I  l.it,  (4.3.4b) 


(4.3.4c) 


and 


I  C*'+i(CT)  I  5 


I.  ae  (-00. 

max  [  I  cU  ( )  I .  i  cU  (  )  I  ]  a  €  ( Xt'^\  Vil  ].  /  =  1 . *-1; 

I  Ci+i(Xt<*^)  I.  06[X*<*\oo). 


Direct  consequences  of  Lemma  4.7  and  Theorem  2.3  in  section  2.4. 


Discussion: 

Lemma  4.5  tells  us  that  c^^i  is  like  the  characteristic  polynomial  of  Tt  in  that  it  vanishes  at  the 
eigenvalues  of  T* .  Moreover  it  is  altematingly  concave  upward  and  concave  downward  in  the  intervals 
divided  by  the  eigenvalues  of  T*.  The  direct  result  from  this  geometric  property  of  Ct+i  is  that 


I  CtVl 


dk^i  ^  U* 


for  a  * 


(4.3.5) 


as  shown  in  the  figure  Fig.  4.2. 
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The  inequality  (4.3.5)  can  be  changed  into: 

H*  I  ct+i  I  S  I  ct+i  I  <  1. 
However,  this  bound  can  be  improved. 


for  any  real  a. 


Lemma  4  For  any  reala, 


p*  I  I  <  1/2, 
p*  I  St^i  I  <  1/2. 


(4.3.6) 

(4.3.7) 


Proof:  When  a  =  then  p^  vanishes  and  the  assertions  are  vacuously  true.  Now  we  suppose  that 


2 

I  '  I  ^*+1  I  '  1 

11*  >  C*+l  I  =  T — p*  1  Jt*  1, 
„S*+1 

<  t —  1  n*  >. 

s*+i 


-  •**+!  I  C*+l 
5  1/2. 


Similarly 


,  -  ,  '  ^t+lC*+l  I  . 

M*  I  5*+i  I  =  - 1 - 11*  I  11*  I . 


^*+1 
S*+l 

=  I  1 , 

<1/2. 


by  (D) 

using  (4.2.9b) 
by(B) 

since  I  s*+iCt+i  I  £  1/2 


by  (H) 
by  (4.2.9b) 


by(B) 

since  I  St+jCt+i  1  S  1/2  and  I  Ct+i  I  <  1. 


As  before  we  need  to  complement  the  proceeding  result  with  another  that  involves  the  eigenvalues  of 

Tt^. 

Lemma  4.9:  For  any  real  a 


14+t  I  <i*+i  I  <  2 


(4.3.8) 


Proof:  When  a  =  then  p*+i  vanishes  and  the  assertion  is  vacuously  true.  Now  we  suppose  that 
a  ^  As  before 


I  it*+i«i*+i  1  =  11+  lt*+lC*+l  I, 

S  1  +  I  Jlt+i  I  I  c*+i  I, 

I  Jt*+i  I  ,  , 

I  Ct+I  I. 


<  I  + 


ltt+1 


by  (E) 
using  (4.2.9b) 


Multiply  each  side  by  Pn+i/ 1  n*+i  I  to  get 


4*+i  I  ^*+1  I  < 


14+1 


Jt*+i  I 
<1+1  c*+i  I  <  2, 


+  I  C|i+i  I, 


by  (G). 


The  above  lemmas  are  summarized  into  the  following  theorem. 
Theorem  4.4:  For  any  real  a 


I  c*+i  I  <  min  [ 


1 


2  p*  ’  P*+i 


]. 


(4.3.9) 
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From  Theorem  4.1  and  Theorem  4J  we  see  that  derivatives  of  at  can  be  huge  if 

I  (Oj^  I  is  tiny  or  I  I  is  tiny  or  I  (Oijc  I  is  tiny. 


4.4  Properties  of  tjk 

This  subsection  concentrates  on  t*  =(Jit(tT) +  Ct (a) pt+i )'.  Since  lit*  II  =  ( ( Jc*  )^  +  ( c*  P*+i 
the  results  in  sections  4.2-4.3  can  be  applied  here  to  bound  1 1 1*  1 1  as  the  following  theorems  show. 

Theorem  4.5:  For  any  real  a 

M^lltill^  <  [|i*  +  -| spread {Tt)f  +  spread\ r*+i ).  (4.4.1) 

Pt^.illtill^  <  j  [  pt. I  +  spread  spread\Tt^i).  (4.4.2) 

Proof;  Use  (4.2.15)  and  (4.3.8)  to  get 

\xi\nt\^  <  [p*  +  ^spreadiTt)]'^  and  4*^lc*'l^  <  4. 

Therefore  apply  Theorem  22  to  get  (4.4.1). 

Use  (4.2.18),  (4.2.10)  and  (4.3.6)  to  get 

IJt*  1^  <  +  spread{Tt)f  and  I  c* 

Apply  Theorem  12  to  get  (4.4.2).  ■ 

Theorem  4.5  shows  that  when  a  is  far  from  the  spectrum  of  T*  and  r*_i,  lit*  II  has  modest  magni¬ 
tude  even  when  the  spectra  of  T*  and  r*_i  have  elfcments  that  are  very  close. 

Remark  1:  A  similar  result  on  II  t*  II  to  Theorem  4.1  and  Theorem  4.3  could  be  derived.  However  it  is  too 
complicated  since  I  n*  I  is  related  to  the  spectra  of  T* ,  r*_i  and  I  c*  I  is  related  to  the  spectra  of 
r*-i,  r*_2  and  the  magnitude  of  P* .  We  present  a  simplified  one. 


Theorem  4.6:  The  norm  of  vector  t*  at  eigenvalues  of  Tt  satisfies 


<  [Jt;(Xi»J)]2  <  iit;(xf))ii, 

<  itti(Xi<*>)]^  <  iit;(x/*))ii,  i  =  2 . k-i 


<  [Jt*(xD]2  <  iu;(xi*>)ii. 


Remark  2:  It  will  be  shown  in  section  5  that  forward  instability  can  appear  at  step  k  only  if  o  is  very  close 
to  an  eigenvalue  of  T* . 

5.  TQR  in  Finite  Precision  Arithmetic 

This  section  studies  the  relations  among  the  computed  quantities  generated  by  TQR  using  finite  pre¬ 
cision  arithmetic.  A  central  objective  is  to  explain  the  different  phenomena  we  have  observed  in  the  exam¬ 
ples  presented  in  section  23.  It  turns  out  that  the  magnitude  of  the  ex^t  II  t*  II  I  sin  ( t* ,  t* )  I  (  defined 
in  Corollary  1  of  Lemma  4.2  )  governs  the  accuracy  of  the  computed  lit*  II,  and  hence  the  accuracy  of 
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the 

computed  nt  and  ct.  From  a  computational  point  of  view,  the  role  of  Hull  i  sin  (U>U)  I  can  be 
replaced  by  the  magnitude  of  the  exact  MU  H  (  to  be  shown  in  (S.2.6) ).  The  combination  of  analysis  in 
this  section  with  Theorem  4,5  tells  us  that  to  have  forward  instability,  it  is  necessary  and  sufficient  that  a 
be  simultaneously  close  to  an  eigenvalue  of  and  an  eigenvalue  of  for  some  k  <n. 

5.1  Perturbed  commutative  law 

We  analyze  TQ/?  in  finite  precision  arithmetic.  If  v  is  an  output  of  TQR  in  exact  arithmetic  then  let  v 
denote  the  corresponding  output  in  finite  precision  arithmetic.  In  particular 


Definition  5.1:  We  denote  by  the  following  vector  in  R* , 

CiC-Ja)  •  •  •  (-U) 
C2(-S3)  •  •  •  (-Sk) 

y*:  = 


(5.1.1) 


Ck-l(-Sk) 

Ck 


It  is  not  the  case  that  y^  and  are  the  ouQ)ut  of  TQR  acting  on  a  perturbed  matrix  with  the  same 
shift  cr.  Nevertheless  y^  and  ttt  do  satisfy  exactly  a  perturbed  version  of  the  fundammtal  relation 

(T*  -  o/*)yk  =  Jttei*>.  /t  =  l . /I.  (5.1.2) 

It  turns  out  that 

( r*  -  OU  +  F*  )y*  =  n* ei*),  k^l . n.  (5.1.3) 

where  Ft  is  a  tridiagonal  matrix  which  is  a  smaU  perturbation  ( component  by  component )  of  T^.  In  fact 

II  Ft  II  <  5.6spreadz,  spread  =  ^(F*)  -  Kiadk)  (5.1.4) 

provided  that 

KJTk)  ^  o  ^  K»^(Tk). 

For  verification  of  relations  (5.1.3)  and  (5.1.4),  see  section  S3. 

The  significance  of  (5.1.3)  lies  in  the  following  perturbed  commutative  law. 

Lemma  5.1:  For  any  real  a, 

c*(a)P*+iJt*(o)  -  Ji*(o)c*(a)P*+,  =  Pt+iy*Fty*.  lS*5n-l.  (5.1.5) 

Proof:  For  simplicity,  we  drop  reference  to  a.  Apply  y*  to  both  sides  of  (5. 1 .3) 

y*(F*  -  <T/ik)y*  +  yt'Fty*  =  ittCt. 

Use  (5.1.2)  and  multiply  by  Pt+i  to  get  (5.1.5).  ■ 

5.2  Forward  instability  of  TQR 

In  section  23,  we  saw  that  TQR  is  forward  unstable  in  Example  2.2  and  forward  stable  in  Example 
2.4.  The  interesting  case  is  Example  2.4.  Since  the  shift  a  =  -2  is  an  eigenvalue  of  Tj  ( the  3x3  leading 
principal  submatrix  of  T5 ),  the  computed  «3(-2)  will  be  tiny.  Obviously,  it  has  lost  all  significant  digits. 
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However  the  computed  }C4(-2)  still  preserves  most  of  its  significant  digits,  unlike  Example  2,2  in  which  the 
computed  KtO^  loses  all  the  significant  digits  for  k  s  4, 5  and  6,  ( see  table  at  the  end  of  section  3.2 ).  The 
key  difference  between  these  two  examples  is  the  lengths  of  residual  vectors  U(-2),  k  =  1, ....  5  in  Exam¬ 
ple  2.4  and  the  lengths  of  residual  vectors  k  =  1, ....  6  in  Example  2.2  as  shown  in  the  following 
table: 


k  lltt(-2)ll 

{Example  2.4) 


{Example  2.2) 


1 

2 

3 

4 

5 

6 


l.OOOOOOOOd+00 

l.OOOOOOOOd-KX) 

5.00000000d-01 

5.000000(X)d-01 

1.30404754d-33 


3.64965820d400 
3.16227602d-03 
9.3588227  Id -07 
2.37408198d-10 
4.48883862d-14 
1.60806628d-17 


The  lengths  of  lltt(-2)ll  never  go  under  the  magnitude  of  0 {spread  ^ )  except  at  last  step.  This  is 
what  we  expected  for  deflation  to  occur.  However  the  length  of  II  tt(X^  II  has  gone  under  the  magnitude  of 
0  {  spread  Ve )  starting  firom  step  4.  The  role  of  II  II  is  shown  by  the  following  relation. 


Corollary  of  Lemma  5.1:  For  any  real  a 

lltt(<i)IIIITt(a)ll  lsin^(tt,Tt)l  =  I  ytT^y*  I,  k  =  l . n-1.  (5.2.1) 

whCTc  T*(a)  =  ( S*.  )' 


Proof:  Use  definitions  of  and  and  observe  that  the  left  hand  side  of  (5.1.5)  is  the  cross-product  of 
these  two  vectors.  Therefore  the  absolute  valiw  of  the  left  hand  side  of  (5.1.5)  is  equal  to 

lltt(<r)ll  IIT*(<t)II  lsin.ti(tt.Tt )  I.  ■ 

A  variation  of  relation  (5.2.1)  is 

llt*(a)ll  lll*(a)-T*(<r)IMsin^(t*.tt-T*)l  =  I  y^F^y*  I.  k  =  l . n-1.  (5.2.2) 


The  combination  of  relation  (5.2.2)  with  Corollary  1  of  Lemma  4.2  that 

lit*  II  lit*  II  I  sin  ^ (t*,  t* )  I  =  P*+i 
yields  the  following  important  relation: 

Theorem  5.1:  For  any  real  a  and  k  =  I, n-l 

llt*(o)-T*(a)ll  lsin..:(t*,t*-T*)l  =  llti II  I  sin  _(t*,ti)  I  I  y^F*y*  I. 

The  direct  consequences  of  Theorem  5.1  are: 

II t* (a) -T* (0)11  2  llt*ll  I  sinz(t*,ti)  1 1  yff*y*  I.  k  =  l . n-1, 

P*+i  I  ylfkfk  I 


and 


lit;  II  s 


lit*  II 


llt*(o)-T*(a)ll  Isin4(t*,t* -T*  )l 


I  VkFkJk  I 


k  =  1 . n-1. 


(5.2.3) 

(5.2.4) 

(5.2.5) 

(5.2.6) 

(5.2.7) 


Now  we  come  back  to  relation  (S.2.1) 

Since  I  P*+i  y*F*y*  I  /2  is  the  AREA  of  the  triangular  defined  by  t* ,  T*  and  t*  -  T* .  We  shall  describe 
instability  in  terms  of  this  geometric  figure.  There  are  three  stages: 


a)  stability 


J  4v, i'l 


b)  onset  of  instability 


characterized  by  relative  error  in  II  II  such  that 
lit* -T*  II  lit* -T*  II 


lit*  11 


IIT*  II 


=  1. 


c)  full  instability 


\1  Hi-  4i  H 


\\4v.\\ 


We  can  expect  to  lose  all  correct  digits  in  the  greater  of  I  a*  I  and  I  c*  I  ^*4.1  when  the  triangle  0 .  t* , 
T*  is  equilateral: 

I  -  I 


^  lit* II"  =  =  p*4iiyi^y*i. 


Example  5.1: 

The  matrix  is  the  same  as  that  in  Example  2.2  with  the  same  shift  The  lengths  of  the  residual  vectors 
from  TQR  at  step  k  are  presented.  For  comparison,  the  correct  lengths  are  also  presented  to  show  the  rela¬ 
tionship  between  the  accuracy  of  IIT*  II  and  the  magniuide  of  lit*  II .  It  is  not  h^  to  observe  that  the  pro¬ 
duct  of  corresponding  elements  in  first  column  and  third  column  is  less  than  spread  e. 


llt*(X«)ll 


IIT*(X«)II 


llt*(X«)-T*(X,i)ll 


1  3.6496582(W+00  3.64965820d-K)0  O.OOOOOOOOd+OO 

2  3.16227602d-03  3.16227602d-03  4.15999586d-17 

3  9.35882271d-07  9.35882271d-07  2.42128268^/-!  3 

4  2.37408198d-10  8.51726993d-10  8.18091259d-10 

5  4.48883862d-l4  3.22498160d-06  3.22498160d-06 

6  1.60806628<i-l7  1. 70563083d -02  1.70563083d -02 
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k 

Isin  1 

dt 

1 

e 

1 

0.0000(W+00 

l.OOOOOfif+OO 

2 

5.26836d-09 

l.OOOOOd+00 

3 

2.58686d-07 

l.OQOOOd-^00 

4 

9.60509d-01 

l.OOOOOd+00 

5 

1.0000(W+00 

l.OOOOOd+00 

6 

0.000(XW400 

O.OOOOCW-fOO 

Example  5.2: 

Another  interesting  example  is  one  sweep  of  TQR  on  the  following  matrix: 

11  11' 

Tjs  =  tridiagonal  13  0  0  •  •  •  0  0  13  . 

11  11 

The  shift  is  its  largest  eigenvalue.  We  exhibit  the  11^11,  lITtll,  I  sin  )  I  and 

I  sin  - 1* )  I  in  the  following  table. 


1 

1.0030d+00 

1.0030d+00 

O.OOOOd+00 

O.OOd+00 

l.OOd+00 

2 

7.6923d-02 

7.6923d-02 

O.OOOOd+00 

0 .OOd+00 

l.OOd+00 

3 

5.9171d-03 

5.9171d-03 

3.0000d-14 

O.OOd+00 

9.97d-01 

4 

4.5516d-04 

4.5516d-04 

3.7620d-13 

0,00d+00 

9.88d-01 

5 

3.5012d-05 

3.5012d-05 

4.8893d-12 

1.38d-07 

9.88d-01 

6 

2.6932d-06 

2.6933d-06 

6.3561d-ll 

2 .33d-05 

9.88d-01 

7 

2.0717d-07 

2.0730d-07 

8.2629d-10 

3.94d-03 

9.88d-01 

8 

1.5936d-08 

2.0536d-08 

1.0742d-08 

5.17d-01 

9.88d-01 

9 

1.2259d-09 

1.3984d-07 

l,3964d-07 

9.87d-01 

9.88d-01 

10 

9.4298d-ll 

1.8154d-06 

1.8154d-06 

9.88d-01 

9.88d-01 

11 

7.2532d-12 

2.3600d-05 

2.3600d-05 

9.88d-01 

9.88d-01 

12 

5.5304d-13 

3.0680d-04 

3.0680d-04 

9.97d-01 

9.97d-01 

13 

5.5304d-13 

3.9883d-03 

3.9883d-03 

7 .67d-02 

7 .67d-02 

14 

7.2532d-12 

5.1848d-02 

5.1848d-02 

4.50d-04 

4.50d-04 

15 

9.4298d-ll 

6.7313d-01 

6.7313d-01 

2 .66d-06 

2.66d-06 

16 

1.2259d-09 

7 .2659d+00 

7.2659d+00 

1.67d-08 

1.58d-08 

17 

1.5936d-08 

1.2916d+01 

1.2916d+01 

O.OOd+00 

O.OOd+00 

18 

2.0717d-07 

1.2999d+01 

1.2999d+01 

0 .OOd+00 

O.OOd+00 

19 

2.6932d-06 

1.3000d+01 

1.3000d+01 

5.27d-09 

5.27d-09 

20 

3.5012d-05 

1.3000d+01 

1.3000d+01 

5.27d-09 

O.OOd+00 

21 

4 .5516d-04 

1.3000d+01 

1.3000d+01 

O.OOd+00 

0 .OOd+00 

22 

5.9171d-03 

1.3000d+01 

1.3006d+01 

O.OOd+00 

0 .OOd+00 

23 

7.6920d-02 

1.3000d+01 

1.3077d+01 

O.OOd+00 

O.OOd+00 

24 

9.9704d-01 

l,3000d+01 

1.3997d+01 

0 .OOd+00 

0 .OOd+00 
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The  behavior  of  II  II  has  gone  through  the  three  stages  described  before  and  the  instability  occurs 
at  step  8.  However  starting  firom  step  13,  the  exact  lift  II  begins  to  increase  and  the  computed  IIT^  II  still 
keeps  increasing.  As  a  result,  the  angle  between  these  two  vectors  has  to  shrink  in  order  to  satisfy  the  rela¬ 
tion  (S.2.1).  Finally  the  triangular  turns  into  the  following  form: 


The  interesting  point  of  this  example  is  that  the  computed  IIT24II  can  be  totally  wrong  even  though 
exact  1 1 124 II  is  not  tiny.  However,  the  instability  occurred  when  I  ITg  1 1  was  about  O  {spread  e). 


Explanation  of  forward  Instability 

The  Corollary  of  Lemma  S.l,  coupled  with  the  behavior  of  ytFgy*,  shows  why  instability  must 
occur  in  certain  cases. 

There  are  classes  of  tridiagonal  matrices  and  shift  a  such  that  II  tg  II  can  diminish  with  il  to  as  small 
a  value  as  desired.  Furthermore  fcv  small  enough  k  the  vectors  and  yt  are  nearly  identical.  As  k 
increases  the  last  components  ct  and  ci  may  decrease  to  0(Ve),  at  which  size  cg  could  lose  all  its 
significant  digits.  Nevotheless  y*  is  an  unit  vector  and  yiyt  =  1  ( i.e.  >  0.8  ).  At  this  stage  ylFiJi  is  close 
to  the  Rayleigh  quotient  of  Ft  for  y^g,  namely  yiP^yt.  For  typical  roundoff  siniations  this  Rayleigh  quotient 
will  be  of  the  same  order  as  IIF^  II .  Recall  ftom  Corollary  of  Lemma  5.2  that  IIFg  II  <  5.6 spread (Tt)e 
when  Kun(Tt)<<s^X^(Ti). 

CONCLUSION: 


a)  If  I  sin.<'(tt,T* -t*)  I  2  V3/2  and  lit*  II  2  4 (Pt+, spread then 

iITt  -  t*  II _ Pt^-i  I  JkPkyk  I  ^  5. spread z  2  ^ 

llt*ll  "  llttll^l  sin^(tt,Tt-t*)  1  163t+, spread e  ^  ^  2’ 


b)  While  I  y*F*y*  I  >  ^  spread  e,  we  have 

111*  -  tk  II  ^  Pt-n  I  yjPkyk  I  ^  P^^■l  spread E 

null  lltgll^l  sinz(t*,Tt-tt)  I  4  111*11^ 


Thus  as  1 1 1*  I P  drops  below  ( P*+i  spread  e  )/4  so  must  the  relative  error  in  T*  rise  drastically  above  1 .  This 
is  complete  instability. 

Note:  we  know  of  one  special  case  ( see  [  Demmel  and  Kahan  ] )  when  all  diagonal  elements  of  T  are  zero 
and  a  variation  of  TQR  with  shift  as  zero  forces  this  condition  on  the  diagonal  elements  of  the  transformed 
T.  In  this  case  it  turns  out  that  ytFgy*  is  structurally  zero  for  all  .  Thus  we  can  not  prove  that  y^F*  y*  in 
general  remains  above  spread  eJ4  until  complete  instability  (  yty*  =  0  ('Jz) )  sets  in.  However  this  is  what 
we  observe. 

5.3  Floating  point  version  of  TQR 

In  the  following,  we  want  to  present  a  floating  point  version  of  the  useful  relation  (3.3.8): 

(T*  -  <J/*  )y*  =  Jt* 

Our  model  of  error  in  floating  point  arithmetic  is 

fl(xoy)  =  (xoyXl+f) 


(5.3.1) 
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wheie  o  is  one  of x  and /;  J?(x  o>  )  is  the  floating  point  result  of  the  operation  o.and  Ic  I  S  e. 
where  e  is  the  machine  precision.  Our  analysis  will  be  linearized  using  Wilkinson's  notation  ( see  [  Wilkin¬ 
son,  pp  113] ): 

If  \  Ci  \  <t,  i  =  l,...,/i  and  ne  S  0.01,  then 

n  ( 1  +  e.- )  =  1  +  1.01  One  (5.3.2) 

i>l 

where  10  1  SI. 


Lemma  5.2;  Let  Ci ,  s,-  and  n,-  denote  the  outputs  of  TQR  in  exact  arithmetic  with  shift  a.  Let  ci ,  Si  and  Ki 
denote  the  outputs  o/TQR  in  finite  precision  arithmetic  with  shift  a.  Then  the  computed  outputs  from 
TQR  with  shift  a  satisfy  exactly  the  following  relation: 

(r*  -  o/*  +  f*)y*  =  S*e*(*)  (5.3.3) 


where 


Ft 


and 


(ai-o)ei,i  pjeu 


Pt6*-ljk 

PtE*jk-i  («t-CT)etjk 


I  e*jk  I  S  3.03e.  I  I  S  3.036,  I  e*jk+i  I  S  2.02e,  for  it  =  1,2 . n. 


Proof:  We  use  mathematical  induction. 

Use  definition  of  tCj  in  (4.1.1)  to  find  rtj  =  ai  -  o.  Hence  by  (5.3.1) 

«i  =  fi(ni)  =  fi(ai-a)  =  (a,-o)(l+e,,i).  (5.3.4) 


Use  definition  of  n*  in  (4.1.1)  and  then  (5.3.1),  (5.3.2)  to  find 
ttt  =  fii  P*Ct-i(-ji)  +  (  ctt  -  CT  )Ct  ) 

=  Pt(  1  +  )ct-\(-st)  +  ( a*  - ct)(  1  +etjt  )ct,  (5.3.5) 

with  I  I  S  3.03e,  I  e*jk  I  S  3.03e  because  each  quantity  is  involved  in  3  atomic  arithmetic 
operations.  Do  the  same  thing  for  ct.st, 

a,-  ,e  ,  ’**-1  -  ara  \  P*(^+^J*) 

Ct  =  =  Y~r, - T-  St  =  fii^kl^k)  =  — = - • 

^*(l+ec*)  %k 

Then 

^kttk—l  ~  P*^*(  1  ^ek  X  1  )  ~  Pt^t(  1  ^k-ljt  )>  (5.3.6) 

with  I  tk-ijc  I  S  2.02e. 

Now  we  start  to  put  things  into  matrix  form.  Consider  the  case  k  =2.  Multiplying  (5.3.4)  by  72, 

jaJti  =  ( Oi  -  a )( 1  +  e,j  )c  i72  since  c  i  =  1 .  (5.3.7) 


Use  (5.3.6)  and  (5.3.5)  for  k  =  2, 
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Jjjjti  =  p2(  1  +  £u_)c2_*  ( Oi  -  o  X 1  +  ei.i  Xif  2; 
«2  =  p2(  1+22.1  )Cl(-J2)  +  (  02-0X1 +E2^)c  2- 
That  is  to  say 


(ai-oXl  +  ei.i)  PXl+eu) 

p2(l  +  e2.i)  (02-0X1  +  22.2) 


by  (5.3.7) 


N 

Ci(-S2) 

0- 

J 

C2 

.  4 

(5.3.8) 


The  fonnulae  (5.3.4)  and  (5.3.8)  have  shown  that  (5.3.3)  is  true  for  ik  =  1, 2.  Now  suppose  the  state¬ 
ment  is  true  for  A:  =  y .  That  is  to  say 

(ai-<T)(l+Ei.,)  PX1+2i2) 

P2(1+22,i)  (0(2-O)(l-f€2j)  P3(l+C23) 

P3(1+23^ 

Py(l+ey_i^) 

Py(l+£yy_,)  (Uj  -OXl+Eyy) 


(5.3.9) 


Ci(-S2)  •  •  *  (-^y) 
C2(rh)  •  •  •  i-Sj) 

f  N 

0 

Vi  J 

0 

.  . 

Multiply  by  (-5}+i)  on  both  sides  of  the  above  equation  and  use  the  results  (5.3.6),  (5.3.5)  fw 
k  =  y+1,  namely 

(— '^>+1)^/'  —  —  P/+l(  1  +  2yy+i  )cy+i> 

P/+l(  1  +  2/+1./  +  (  Oty+2  —  0)(  1  +  Ey+iy+l  )c^+l  *  tty.,.1. 


Now  adjust  the  last  row  of  (5.3.9)  and  add  a  new  row  to  get 


(ai-CT)(l+ei,i)  P2(l+eij) 

P2(l+^2.l)  (02-<T)(  1+62.2)  P3(1+224) 

P3(  1+23.2) 


P>(l+2y-l.y) 

P/(l+e7,/-l)  (a^-cXl+Eyy)  P/+l(l+eyj>i) 

P/+l(l+ey+,^)  (ay+i-C)(l+€y+,^>,) 


* 

Ci(-J2)  ■  •  •  (-Jy+i) 

f  > 

0 

C2(-S3)  •  •  •  (-i>+i) 

^y(— ^/+i) 

0 

Cy+i 

•i  t 

,  4 

By  induction,  (5.3.3)  holds  for  it  ^  n .  ■ 
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CoroUary  of  Lemma  5.2  If 


then 


X^iT^)  ^  a  ^  KJiT,) 


II  F,  II  <  5.6  spread 

where  spread T,  )  -  X^( T,  ). 


(5.3.10) 


Proof: 

(1)  by  Lemma  5.1 

I  etjk  I  5  3.03e,  I  c*jk_i  I  5  3.03e.  I  etjk+i  I  ^  2.02e,  for  k  =  1,2 . n; 

(2)  I  I  S  spread{T^)l2  for  2^k^n,  by  Theorem  2J!; 

(3)  I  ct*  -  <J  I  5  spreadiT^)  for  by  condition  given, 

then,  abbreviating  spread  ( T,  )  by  spread , 

II  F,  IL  =  itm  [  I  pt  I  +  I  ( ot*  -  a)  etjk  I  +  I  Pt+,  e^jk+i  I  ] 

3  2 

S  —  spread  1.01e  +  3  spread  l.Ole  +  —  spread  l.Ole 

=  (3/2 +  3  +  1)  1.01  spread  t 

<  5.6  spread  e. 

In  the  same  way, 

II  F,  llj  <  5.6  spread  e. 

Hence 

II  F.  IP  =  inaxX((F,)'F,) 

5  II  (F,)‘F,  IL 
5  II  (F,  )'  II.  II  F,  II. 

=  II  F,  II,  II  F,  II. 

<  ( 5.6  spread  e  )*.  ■ 

Note:  When  <T  is  outside  the  interval  fX..nfr.). X^-.fT.)!  we  can  bound  term  I  a* -a  I  by 
ji,  +  spread  (T*).  Therefore  for  any  real  a 

IIF,  II  <  5.6spread(TJt  +  3.03Hne.  (5.3.11) 


5.4  Forward  stability 

In  this  subsection,  we  will  present  a  forward  perturbation  analysis.  It  turns  out  that  the  instability  can 
occur  only  if  a  is  very  close  to  an  eigenvalue  of  . 

Definition  5.2: 

J*:  =  rtt:  =  jtt/lly*ll. 

Recall  that  p»(a)  =  mini  <t-  A,/*’  I  and  spread{T^)  =  X^(Tt)-  X^{Tt). 

i 

Theorem  5.2  For  any  real  o,  if 

Pt(a)  >  12j/;rea<f(ri)£,  t  =  1 . n  (5.4.1) 


then 
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^ky  - 


I  wt  -  I  ]_ 

I  n*  I  ^  2  ’ 


Proof:  Let  B*  =  (r*  -  a/*  r'F*.  Use  (5.1.3)  to  find 


k  =  l . n. 


Sk  -  (.Tic-alh+Fh)  ‘etrtt. 


(5.4.2) 


=  (/+B*)-‘(r*-o/*)-‘e*7i*. 

=  (/+Bt)-‘y*nt/Jtt.  (5.4.3) 


Premultiply  by  yi  to  get  the  main  formula 

Jkh  =  y*(/+Bt)~‘y*n*/Jt*. 

(1)  we  show 

IIB*II  S  1/2,  ifc  =  l,2.....n. 
IIF*1I 


IIB*II  5  ll(Tt-a/ir  II  ll^tll  S 
3.03  p*e  +  5.65/»refld(r*)£ 


14 


<  1/2, 


14 


by  (5.3.11) 
by  (5.4.1).  ■ 


(2)  we  show 


yi(/  +  Si  )"'yt  >0.  i  =  1, 2 . n. 


(5.4.4) 

(5.4.5) 


(5.4.6) 


y*'(/+F*)-'y*  =  1  +  i(-i)'y*'5iy* 

1*1 

s  1  -  £  iiBiii, 

i  9  1 

>  1  _  - 5 - 

1  -  IIB*II  ’ 

1  -  2  IIBtll 

"  1  -  IIBtll  ’ 

>  0,  by  (5.4.5). 

(3)  we  show  by  induction  that 

y*y*  >  0.  Titn*  >  0,  *  =  1,2 . n.  (5.4.7) 

Note  that  y{ft  =  1  >  0  and  WiJii  >  0  by  (5.4.4)  and  (5.4.6)  with  *  =  1.  Now  suppose  that 

yjS/  >  0,  n^Kj  >  0,  (5.4.8) 


then 


^/+l  4y+l 


(1  +£)  = 


Wy  II  y>  II 
^y+i  ^y+i 


(1+e)  >  0. 


Reference  to  (5.1.1)  shows  that 


(5.4.9) 


yy+i^y+i  =  yy+iyy+i/iiyy+iii  =  ".r:— "  (yjyy  “yy Wy>i  +  c^+icy+i). 

Ilyy+ill 

Since  sj+t  >  0,  sj+i  >  0  ( because  Py+i  >  0  )  then  by  (5.4.8)  and  (5.4.9)  we  see  that 

yy+iJy+i  >  0. 

Use  (5,4.4)  and  (5.4.6)  with  k  =  j+\  to  find  Jty+iny+i  >  0.  By  principle  of  induction,  (5.4.7)  is  uiie 


-35- 


fori  =  1, .... /I. 

(4)  we  show  (5.4^). 

Take  norms  of  (5.4.3)  and  use  (5.4.7)  to  get 

Ttt/jr*  =  II  (/  +  5*)"‘y*  <1 

Use  standard  perturbation  theory  and  (5.4.5)  to  find: 

Ii(/+5tr‘y*  II  ^  Ii(/+fit)-‘ii  s 


Let  y*  =  U  +^k )w.  Then  1^(1+  MB*  II )  llwll.  Hence 

- ^ 5  II  (/ +Bi)-'yt  II. 

1+  IIB*I1  ^ 


Therefore 


So 


1  -  1 
<  < 


1  +  IIBk  II 


1-  IIBtll 


I  «*  -  «*  I  ^  »  M  1 

=  Irt*  I  ^  ^  2- 


Conclusion 

Forward  instability  was  introduced  as  the  occurrence  of  a  computed  transform  f  that  was  far  from 
the  exact  f.  However  our  study  concentrates  attention  on  the  forward  stability  of  the  rotation  angles  that 
determine  the  similarity  transftsrmation.  This  is  reasonable  because  wildly  oroneous  angles  ensures  both  a 
wrong  f  and  also  a  vector  y,  that  is  not  a  reasonable  eigenvector  approximation.  However  even  rotation 
angles  with  high  relative  accuracy  cannot  guarantee  that  each  computed  element  of  7  has  high  relative 
accuracy. 

In  section  5.4  we  have  shown  that  the  computed  will  not  ^ose  all  the  significant  digits  if  the  shift  a 
is  nowhere  close  to  any  eigenvalues  of  Tj,  j  =  \,  That  is  equivalent  to  say  that  instability  can  only 
happen  at  step  k  if  a  has  been  very  close  to  some  eigenvalue  of  Tj,j  <k.  However  a  being  close  to  an 
eigenvalue  itself  does  not  provoke  forward  instability.  Example  2.4  illustrated  this  fact. 

We  point  out  some  main  results  of  this  study. 

(1)  The  effect  on  and  of  changes  in  the  matrix  elements  Oi  and  is  transmitted  through 
a  perturbation  matrix  .  By  Theorem  5.1 

IIT*-t*lllsini.(tt.Tt-U)l  <  null  lyiF*y*l. 

and  the  amplification  factor  for  the  lytFjyt  I  is  lltj  II,  the  derivative  with  respect  to  o. 
Instability  requires  that  11 U  H  be  huge.  There  is  no  need  to  compute  the  gradient  of  II U  H 
with  respect  to  all  the  variables  oc,-  and  Pi . 

(2)  Formula  (5.2.6)  shows  that 

llU-Ttll  \ylF,y,  I 

111*11  lit*  11^ 

Thus  whenevo-  lit*  II  gets  too  small,  the  relative  error  in  lit*  II  can  rise  well  above  1.  In 
other  words  small  enough  values  of  Ic*  I  and  lie*  I  cannot  be  calculated  accurately  in  finite 


L 
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precision  arithmetic.  Actually  when  lit*  11  is  tiny,  n*  and  ct^Pt+i  will  be  tiny.  Therefore 
current  shift  o  has  to  be  very  close  to  an  eigenvalue  of  T*.  Then  by  (4.3.3)  and  (4.1.6)  we 
see  that  1/1  CtPt+i  I  will  be  very  close  to  the  derivative  of  c*+i  at  this  a.  Therefore  tiny 
lltt  II  signals  huge  derivative  on  c^+i.  Using  Theorem  4.4,  we  see  that  p*  and  Pt+i  have  to 
be  tiny  simultaneously. 
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QR  is  the  standard  method  for  finding  all  the  eigenvalues  of  a  symmetric  tridiagonal  matrix.  It  produces  a 
sequence  of  similar  tridiagonals.  It  is  well  known  that  the  QR  transformation  from  T  to  f  is  backward 
stable.  That  means  that  the  computed  f  is  exactly  orthogonally  similar  to  a  matrix  close  toT.  It  is  also 
known  that  the  algorithm  sometimes  exhibits  forward  instability.  That  means  that  the  computed  t  is  not 
close  to  the  exact  f. 

For  the  purpose  of  computing  eigenvalues  the  property  of  backward  stability  is  alt  that  one  requires.  How¬ 
ever  the  QR  trandbrmation  has  other  uses  and  there  forward  stability  is  wanted. 

This  report  analyzes  the  forward  instability  and  shows  that  it  occurs  only  when  the  shift  causes  premature 
deflation.  We  show  that  forward  stability  is  governed  by  the  behavior,  in  exact  arithmetic,  of  a  pair  cf  vari¬ 
ables  and  we  establish  tight  upper  and  lower  bounds  on  their  derivatives  with  respect  to  change  in  the  shift 
parameter. 
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